Fast Algorithms for Multidimensional Signals

Similar to 1-D Digital signal processing in case of the Multidimensional signal processing we have Efficient algorithms. The efficiency of an Algorithm can be evaluated by the amount of computational resources it takes to compute output or the quantity of interest. In this page, two of the very efficient algorithms for multidimensional signals are explained. For the sake of simplicity and description it is explained for 2-D Signals. However, same theory holds good for M-D signals. The exact computational savings for each algorithm is also mentioned.

Motivation and applications

In the case of digital systems, a mathematical expressions can be used to describe the input-output relationship and an algorithm can be used to implement this relationship. Similarly, algorithms can be developed to implement different transforms such as Digital filter, Fourier transform, Histogram, Image Enhancements, etc. Direct implementation of these input-output relationships and transforms is not necessarily the most efficient way to implement those.

When people began to compute such outputs from input through direct implementation, they began to look for more efficient ways. This wiki page aims at showcasing such efficient and fast algorithms for multidimensional signals and systems. A multidimensional (M-D) signal can be modeled as a function of M independent variables, where M is greater than or equal to 2. These signals may be categorized as continuous, discrete, or mixed. A continuous signal can be modeled as a function of independent variables which range over a continuum of values, example – an audio wave travelling in space, 3-D space waves measured at different times. A discrete signal, on the other hand, can be modeled as a function defined only on a set of points, such as the set of integers. An image is the simplest example of a 2-D discrete domain signal that is spatial in nature.

In the context of Fast Algorithms, consider the example below:

We need to compute A which is given by

A = αγ + αδ + βγ + βδ where α,β,γ and δ are complex variables.

To compute A, we need 4 complex multiplications and 3 complex additions. The above equation can be written in its simplified form as

A = (α + β)(γ + δ)

This form requires only 1 complex multiplication and 2 complex additions.

Thus the second way of computing A is much more efficient and fast compared to the first method of computing A. This is the motivation for the evolution of the fast algorithms in the digital signal processing Field. Consequently, many of the real-world applications make use of these efficient Algorithms for fast computations.

Problem Statement and Basics

The simplest form of representing a Linear Shift Invariant system(LSI) is through its Impulse response. The output of such LSI discrete domain system is given by the convolution of its input signal and system's impulse response. This is mathematically represented as follows:

[math]\displaystyle{ y\left(n_1,n_2\right)=\sum_{l_1=-\infty}^{+\infty}\sum_{l_2=\infty}^{+\infty} x(l_1,l_2) h(n_1-l_1,n_2-l_2) }[/math]

where [math]\displaystyle{ h(n_1,n_2) }[/math] is the impulse response of the system.

According to the equation above to obtain the Output value at a particular point ( say [math]\displaystyle{ y(0,0) }[/math]) we need to multiply several values of the input [math]\displaystyle{ x(l_1,l_2) }[/math] and Impulse Response [math]\displaystyle{ h(-l_1,-l_2) }[/math]. Of course this is dependent on the Region of Support of the input as well as the impulse response. The key point here to be noted is that we need to perform so many complex multiplications and additions to obtain 1 output value.

Assuming a 2-D input signal is of length [math]\displaystyle{ M \times M }[/math] and the system's impulse response is of length [math]\displaystyle{ N \times N }[/math] we need to perform [math]\displaystyle{ M^2N^2 }[/math] number of multiplications to obtain all output values. The output can be computed efficiently if one can exploit some characteristics of the system.

We encounter a similar scenario when we have to compute the discrete Fourier Transforms of a signal of interest.

The Direct Calculation of the 2-D DFT is simply the evaluation of the double Sum

[math]\displaystyle{ X\left(k_1,k_2\right)=\sum_{n_1=0}^{N_1-1}\sum_{n_2=0}^{N_2-1} x(n_1,n_2) W_{N_1}^{n_1k_1}W_{N_2}^{n_2k_2} }[/math] [math]\displaystyle{ \text{ for } 0\leq k_1\leq N_1-1 \text{ and } 0\leq k_2\leq N_2-1 }[/math]

[math]\displaystyle{ \text{ where } W_N \equiv \text{exp}\left(-j\frac{2\pi}{N}\right) }[/math]

The total number of complex multiplications and complex additions needed to evaluate this 2-D DFT by direct calculation is [math]\displaystyle{ N_1^2N_2^2 }[/math]. This is a naive approach, however, we already know that an N-point 1-D DFT can be computed with far fewer than [math]\displaystyle{ N^2 }[/math] multiplications by using the Fast Fourier Transform (FFT) algorithm. As described in the next section we can develop Fast Fourier transforms for calculating 2-D or higher dimensional DFTs as well

Fast Algorithms for Multidimensional signals

Row Column Decomposition approach for the evaluation of DFT

Source:

The DFT sum [math]\displaystyle{ X\left(k_1,k_2\right) }[/math] in the previous equation can also be written in the following form

[math]\displaystyle{ X\left(k_1,k_2\right)=\sum_{n_1=0}^{N_1-1}\left[\sum_{n_2=0}^{N_2-1} x(n_1,n_2) W_{N_2}^{n_2k_2}\right]W_{N_1}^{n_1k_1} }[/math]

Let [math]\displaystyle{ G\left(n_1,k_2\right) }[/math] denote the quantity inside the brackets and is given by:

[math]\displaystyle{ G\left(n_1,k_2\right)=\sum_{n_2=0}^{N_2-1} x(n_1,n_2) W_{N_2}^{n_2k_2} }[/math]

[math]\displaystyle{ X\left(k_1,k_2\right)=\sum_{n_1=0}^{N_1-1}G\left(n_1,k_2\right) W_{N_1}^{n_1k_2} }[/math]

Employing this method the DFT [math]\displaystyle{ X }[/math] can be computed as multiple 1-D DFTs. That is, each column of [math]\displaystyle{ G }[/math] can be considered as a 1-D DFT of the corresponding column of [math]\displaystyle{ x }[/math]([math]\displaystyle{ n_1 }[/math] = constant). And each row of [math]\displaystyle{ X }[/math] is the 1-DFT of the corresponding row of the [math]\displaystyle{ G }[/math]([math]\displaystyle{ n_2 }[/math] = constant). Hence we are computing the 2-D DFT by decomposing it into Row and Column DFTs.

The same principle is employed for evaluating the M-D DFT of a M - dimensional signal.

Now let's talk about the computational savings we get using this approach. It is observed that we require [math]\displaystyle{ N_1N_2 (N_1 + N_2) }[/math] complex additions and multiplications. Further, if each of these 1-D DFT is computed using a 1-D FFT, the number of complex multiplications can be further reduced to [math]\displaystyle{ N_1N_2\frac{\log_2 N_1N_2}{2} }[/math]

Vector Radix Fast Fourier Transform

Source:

Just like the 1-D FFT, decimation in time can be achieved in the case of 2-D Signals. The 1-D DFT of a signal whose length is a power of 2, can be expressed in terms of two half-length DFTs, each of these can again be expressed as a combination of quarter length DFTs and so on.

In the case of 2-D signals we can express the [math]\displaystyle{ (N_1 \text{ x } N_2) }[/math] DFT in terms of four [math]\displaystyle{ \frac{N_1}{2} \text{ x }\frac{N_2}{2} }[/math] DFTs (assuming [math]\displaystyle{ N_1 }[/math] and [math]\displaystyle{ N_2 }[/math] are powers of 2). For the sake of simplicity let us assume that [math]\displaystyle{ N_1 = N_2 = N }[/math]. The DFT double sum can be decomposed into four separate summations, one over those samples of [math]\displaystyle{ x }[/math] for which both [math]\displaystyle{ n_1 }[/math] and [math]\displaystyle{ n_2 }[/math] are even, one for which [math]\displaystyle{ n_1 }[/math] is even and [math]\displaystyle{ n_2 }[/math] is odd, one for which [math]\displaystyle{ n_1 }[/math] is odd and [math]\displaystyle{ n_2 }[/math] is even and the last one for which [math]\displaystyle{ n_1 }[/math] and [math]\displaystyle{ n_2 }[/math] are odd.

This is written as :

[math]\displaystyle{ X\left(k_1,k_2\right)=S_{00}(k_1,k_2)+S_{01}(k_1,k_2)W_{N}^{k_2}+S_{10}(k_1,k_2)W_{N}^{k_1}+S_{11}(k_1,k_2)W_{N}^{k_1+k_2} }[/math]

where

[math]\displaystyle{ S_{00}(k_1,k_2) \equiv \sum_{m_1=0}^{N/2-1}\sum_{m_2=0}^{N/2-1} x(2m_1,2m_2) W_{N}^{2m_1k_1+2m_2k_2} }[/math]

[math]\displaystyle{ S_{01}(k_1,k_2) \equiv \sum_{m_1=0}^{N/2-1}\sum_{m_2=0}^{N/2-1} x(2m_1,2m_2+1) W_{N}^{2m_1k_1+2m_2k_2} }[/math]

[math]\displaystyle{ S_{10}(k_1,k_2) \equiv \sum_{m_1=0}^{N/2-1}\sum_{m_2=0}^{N/2-1} x(2m_1+1,2m_2) W_{N}^{2m_1k_1+2m_2k_2} }[/math]

[math]\displaystyle{ S_{11}(k_1,k_2) \equiv \sum_{m_1=0}^{N/2-1}\sum_{m_2=0}^{N/2-1} x(2m_1+1,2m_2+1) W_{N}^{2m_1k_1+2m_2k_2} }[/math]

All the arrays [math]\displaystyle{ S_{00} S_{01} S_{10} }[/math] and [math]\displaystyle{ S_{11} }[/math] are each periodic in [math]\displaystyle{ (k_1,k_2) }[/math] with horizontal and vertical periods [math]\displaystyle{ \frac{N}{2} }[/math]. Using this fact and also the fact that [math]\displaystyle{ W_{N}^{N/2} = -1, }[/math] we can obtain the following identities :

[math]\displaystyle{ X\left(k_1,k_2\right)=S_{00}(k_1,k_2)+S_{01}(k_1,k_2)W_{N}^{k_2}+S_{10}(k_1,k_2)W_{N}^{k_1}+S_{11}(k_1,k_2)W_{N}^{k_1+k_2} }[/math]

[math]\displaystyle{ X\left(k_1+\frac{N}{2},k_2\right)=S_{00}(k_1,k_2)+S_{01}(k_1,k_2)W_{N}^{k_2}-S_{10}(k_1,k_2)W_{N}^{k_1}-S_{11}(k_1,k_2)W_{N}^{k_1+k_2} }[/math]

[math]\displaystyle{ X\left(k_1,k_2+\frac{N}{2}\right)=S_{00}(k_1,k_2)-S_{01}(k_1,k_2)W_{N}^{k_2}+S_{10}(k_1,k_2)W_{N}^{k_1}-S_{11}(k_1,k_2)W_{N}^{k_1+k_2} }[/math]

[math]\displaystyle{ X\left(k_1+\frac{N}{2},k_2+\frac{N}{2}\right)=S_{00}(k_1,k_2)-S_{01}(k_1,k_2)W_{N}^{k_2}-S_{10}(k_1,k_2)W_{N}^{k_1}+S_{11}(k_1,k_2)W_{N}^{k_1+k_2} }[/math]

The above equation tell us how to compute the four DFT points [math]\displaystyle{ X\left(k_1,k_2\right) X\left(k_1+\frac{N}{2},k_2\right) X\left(k_1,k_2+\frac{N}{2}\right) \text{ and } X\left(k_1+\frac{N}{2},k_2+\frac{N}{2}\right) }[/math] for a particular value of [math]\displaystyle{ (k_1,k_2) }[/math] from the four points [math]\displaystyle{ S_{00}(k_1,k_2) S_{01}(k_1,k_2) S_{10}(k_1,k_2) \text{ and } S_{11}(k_1,k_2) }[/math]. [math]\displaystyle{ S_{00}(k_1,k_2) }[/math] can be obtained by evaluating a [math]\displaystyle{ \left(\frac{N}{2} \times \frac{N}{2} \right) }[/math] -point DFT (similarly other [math]\displaystyle{ S_{ij} }[/math] can be obtained).

Thus we see that the [math]\displaystyle{ N \times N }[/math] DFT can be expressed in terms of four [math]\displaystyle{ \frac{N}{2} \times \frac{N}{2} }[/math] DFTs.

By analogy from the 1-D case, the computation depicted in the below figure is called a [math]\displaystyle{ \text{butterfly} }[/math] or more precisely [math]\displaystyle{ radix - (2 \times 2)\text{ butterfly} }[/math] .

Each butterfly requires three complex multiplications and eight complex additions for the calculation of outputs from the inputs. And to compute all the samples of [math]\displaystyle{ X }[/math] from [math]\displaystyle{ S_{00}(k_1,k_2), S_{01}(k_1,k_2), S_{10}(k_1,k_2) \text{ and } S_{11}(k_1,k_2) }[/math] it requires calculations of [math]\displaystyle{ \frac{N^2}{4} }[/math] butterflies.

This decimation procedure is performed [math]\displaystyle{ \log_2 N }[/math] times when [math]\displaystyle{ N }[/math] is a power of 2. Each stage of decimation consists of [math]\displaystyle{ \frac{N^2}{4} }[/math] butterflies, and each butterfly involves three complex multiplications and eight complex additions and hence the number of complex multiplications that needs to be performed during the computation of an [math]\displaystyle{ (N \times N) }[/math] -point radix [math]\displaystyle{ ( 2 \times 2) }[/math] FFT is given by

[math]\displaystyle{ C_{vectorRadix(2\times2)} = \frac{3N^2}{4}\log_2 N }[/math]

See also

• Fast Fourier transform
• Multidimensional transform
• Multidimensional sampling
• Multidimensional discrete convolution
• Multidimensional filter design and Implementation